Modelos - Punto Extra

1 Intrucciones

Considera dos variables \(X\) y \(Y\) con distribución exponencial donde una tiene valor esperado que sea el doble del valor esperado de la otra (tú decides el parámetro).

1. Escoge al menos 5 tipos de cópulas para modelar la dependencia entre ellas, para cada cópula escoge parámetros que reflejen una dependencia negativa y una dependencia positiva (en total tendrás 10 modelos). Para cada uno de los 10 modelos:
   1. Muestra la función de distribución conjunta de \(X\) y \(Y\) (la fórmula).
   2. Grafica la función de densidad conjunta de \(X\) y \(Y\).
   3. Varianza de \(X\), \(Y\) y \(X + Y\).
   4. Desviación estándar de \(X\), \(Y\) y \(X + Y\).
   5. Cuantil al 95% de \(X\), \(Y\) y \(X + Y\).
   6. Valor esperado condicional de la cola al 95% de \(X\), \(Y\) y \(X + Y\).
2. Repite el inciso anterior; pero ahora considera que \(X\) y \(Y\) tienen distribución pareto donde una tiene valor esperado que sea el doble del valor esperado de la otra (tú decides los parámetros).

2 Solución

2.1 Caso Exponencial

Tenemos las siguientes variables aleatorias:

\[X\sim \textbf{Exp}(\lambda)\]

\[Y\sim \textbf{Exp}(2\lambda)\]

donde \(F_{x}(x)=1-e^{-\lambda x}\) y \(F_{y}(y)=1-e^{-2\lambda y}\)

Código

library(copula)
library(lattice)

# Exponential distribution for both marginals
lambda <- 1
# Parameters for the marginals
# You can adjust these parameters as needed
params1 <- lambda # Rate for the first exponential distribution
params2 <- 2*lambda # Rate for the second exponential distribution

2.1.1 Cópula Gaussiana

\[C(u,v)=\Phi_{2}(\Phi^{-1}(u),\Phi^{-1}(v), \rho)\]

donde \(\Phi_{2}\) es la función de distribución conjunta de una distribución normal bivariada con media cero y varianza uno y correlación \(\rho\) (esta determinará la relación negativa o positiva de las variables); y \(\Phi\) es la función de distribución de una distribución normal estándar.

2.1.1.1 Función de distribución conjunta

\[H(x, y) = \Phi_{2}\left(\Phi^{-1}(1 - e^{-\lambda x}), \Phi^{-1}(1 - e^{-2\lambda y}); \rho \right)\]

2.1.1.2 Función de densidad conjunta

Código

library(copula)

rho =.9

# Definir cópula gaussiana
myCopulaP <- normalCopula(rho, dim = 2)
myCopulaN <- normalCopula(-rho, dim = 2)

# Número de muestras
n <- 10000

# Generar muestras de la cópula
set.seed(123)
samples <- rCopula(n, myCopulaP) # caso positivo
samples2 <- rCopula(n, myCopulaN) # caso negativo

# Convertir las muestras a distribuciones exponenciales

## caso positivo
samples[,1] <- qexp(samples[,1], rate = params1)
samples[,2] <- qexp(samples[,2], rate = params2)

## caso negativo
samples2[,1] <- qexp(samples2[,1], rate = params1)
samples2[,2] <- qexp(samples2[,2], rate = params2)


library(ks)

# Estimación de la densidad
dens <- kde(as.matrix(samples)) # caso positivo
dens2 <- kde(as.matrix(samples2)) # caso negativo

library(plotly)

# Convertir la estimación de densidad a un formato adecuado

## caso positivo
x <- dens$eval.points[[1]] 
y <- dens$eval.points[[2]]
z <- matrix(dens$estimate, nrow = length(x), ncol = length(y))
## caso negativo
x2 <- dens2$eval.points[[1]] 
y2 <- dens2$eval.points[[2]]
z2 <- matrix(dens2$estimate, nrow = length(x2), ncol = length(y2))

2.1.1.2.1 Caso Positivo

Código

# Crear gráfico 3D
## caso positivo
plot_ly(x = ~x, y = ~y, z = ~z, type = "surface")

2.1.1.2.2 Caso Negativo

Código

plot_ly(x= ~x2, y = ~y2, z = ~z2, type = "surface")

2.1.1.3 Varianza

Código

# Calcula las estadísticas para las muestras generadas
# Para la correlación negativa (o positiva, según corresponda)
## Caso positivo
X_Positivo <- samples[,1]
Y_Positivo <- samples[,2]
Sum_XY_Positivo <- X_Positivo + Y_Positivo

## Caso negativo
X_Negative <- samples2[,1]
Y_Negative <- samples2[,2]
Sum_XY_Negative <- X_Negative + Y_Negative

# Varianza
var_X <- 1
var_Y <- 1/4
## Caso Positivo
var_Sum_XY_Positivo <- var(Sum_XY_Positivo)
## Caso Negativo
var_Sum_XY_Negative <- var(Sum_XY_Negative)

# Desviación estándar
sd_X <- sqrt(var_X)
sd_Y <- sqrt(var_Y)
## Caso Positivo
sd_Sum_XY_Positivo <- sd(Sum_XY_Positivo)
## Caso Negativo
sd_Sum_XY_Negative <- sd(Sum_XY_Negative)

# Cuantil al 95%
quantile_X_95 <- qexp(0.95, rate = params1)
quantile_Y_95 <- qexp(0.95, rate = params2)
## Caso positivo
quantile_Sum_XY_95_Positivo <- quantile(Sum_XY_Positivo, 0.95)
## Caso Negativo
quantile_Sum_XY_95_Negative <- quantile(Sum_XY_Negative, 0.95)

# Valor esperado condicional de la cola al 95%
## Caso positivo
#exp_X_Tail_95 <-
#exp_Y_Tail_95 <- 
mean_Sum_XY_Tail_95_Positivo <- mean(Sum_XY_Positivo[Sum_XY_Positivo > quantile_Sum_XY_95_Positivo])
## Caso Negativo
mean_Sum_XY_Tail_95_Negative <- mean(Sum_XY_Negative[Sum_XY_Negative > quantile_Sum_XY_95_Negative])

# Imprimir los resultados
cat("Correlación Negativa:\n")

Correlación Negativa:

Código

cat("Varianzas: X =", var_X, ", Y =", var_Y, ", X + Y =", var_Sum_XY_Negative, "\n")

Varianzas: X = 1 , Y = 0.25 , X + Y = 0.634196 

Código

cat("Desviaciones Estándar: X =", sd_X, ", Y =", sd_Y, ", X + Y =", sd_Sum_XY_Negative, "\n")

Desviaciones Estándar: X = 1 , Y = 0.5 , X + Y = 0.7963643 

Código

cat("Cuantiles al 95%: X =", quantile_X_95, ", Y =", quantile_Y_95, ", X + Y =", quantile_Sum_XY_95_Negative, "\n")

Cuantiles al 95%: X = 2.995732 , Y = 1.497866 , X + Y = 3.000276 

Código

cat("Valor esperado condicional de la cola al 95%: X =", 2, ", Y =", 2, ", X + Y =", mean_Sum_XY_Tail_95_Negative, "\n")

Valor esperado condicional de la cola al 95%: X = 2 , Y = 2 , X + Y = 3.993384 

Código

# Imprimir los resultados
cat("Correlación Positiva:\n")

Correlación Positiva:

Código

cat("Varianzas: X =", var_X, ", Y =", var_Y, ", X + Y =", var_Sum_XY_Positivo, "\n")

Varianzas: X = 1 , Y = 0.25 , X + Y = 2.139616 

Código

cat("Desviaciones Estándar: X =", sd_X, ", Y =", sd_Y, ", X + Y =", sd_Sum_XY_Positivo, "\n")

Desviaciones Estándar: X = 1 , Y = 0.5 , X + Y = 1.462743 

Código

cat("Cuantiles al 95%: X =", quantile_X_95, ", Y =", quantile_Y_95, ", X + Y =", quantile_Sum_XY_95_Positivo, "\n")

Cuantiles al 95%: X = 2.995732 , Y = 1.497866 , X + Y = 4.390258 

Código

cat("Valor esperado condicional de la cola al 95%: X =", 2, ", Y =", 2, ", X + Y =", mean_Sum_XY_Tail_95_Positivo, "\n")

Valor esperado condicional de la cola al 95%: X = 2 , Y = 2 , X + Y = 5.877844 

2.1.2 Cópula t-Student

\[C(u,v)=\textbf{t}_{\nu}(\textbf{t}_{\nu}^{-1}(u),\textbf{t}_{\nu}^{-1}(v); \rho)\]

donde \(\textbf{t}_{\nu}(;\rho)\) es la función de distribución conjunta de una distribución t-Student bivariada estandarizada y con correlación \(\rho\) (esta determinará la relación negativa o positiva de las variables); y \(\textbf{t}_{\nu}\) es la función de distribución de una distribución univaria t-Student.

2.1.2.1 Función de distribución conjunta

\[H(x, y) = \textbf{t}_{\nu}(\textbf{t}_{\nu}^{-1}(1 - e^{-\lambda x}), \textbf{t}_{\nu}^{-1}(1 - e^{-2\lambda y}); \rho)\]

2.1.2.2 Función de densidad conjunta

Código

source("codigo.R")
# Parámetros de la cópula y las distribuciones marginales
rho = 0.9
df = 4 # Grados de libertad para la cópula t de Student
params1 = 1 # Parámetro para la primera distribución exponencial
params2 = 1 # Parámetro para la segunda distribución exponencial

# Definir cópulas t de Student
myCopulaP <- tCopula(rho, dim = 2, df = df)
myCopulaN <- tCopula(-rho, dim = 2, df = df)

grafCopula2 <- graficaDensidadCopulas(myCopulaP,myCopulaN)

2.1.2.2.1 Caso Positivo

Código

# Mostrar los gráficos
grafCopula2[1]

[[1]]

2.1.2.2.2 Caso Negativo

Código

# Mostrar los gráficos
grafCopula2[2]

[[1]]

2.1.3 Varianza

2.1.4 Arquimediana Clayton

\[C(u,v)=\phi^{-1}(\phi(u)+\phi(v))\]

donde \(\phi(u)=\frac{1}{\theta}\left(t^{-\theta}-1\right)\)

2.1.4.1 Función de distribución conjunta

\[H(x, y) = \max \left\{ \left((1 - e^{-\lambda x})^{-\theta} + (1 - e^{-2\lambda y})^{-\theta} - 1\right)^{-1/\theta}, 0 \right\}\]

2.1.4.2 Función de densidad conjunta

Código

# Definir parámetros de la cópula de Clayton
theta1 = 2 #  dependencia positiva
theta2 = -.9 #  dependencia negativa

# Definir cópulas de Clayton
myCopula1 <- claytonCopula(theta1, dim = 2)
myCopula2 <- claytonCopula(theta2, dim = 2)

# Número de muestras
n <- 10000

# Generar muestras de la cópula
set.seed(123)
samples1 <- rCopula(n, myCopula1) #  dependencia positiva
samples2 <- rCopula(n, myCopula2) #  dependencia negativa

# Convertir las muestras a distribuciones exponenciales
samples1[,1] <- qexp(samples1[,1], rate = params1)
samples1[,2] <- qexp(samples1[,2], rate = params2)
samples2[,1] <- qexp(samples2[,1], rate = params1)
samples2[,2] <- qexp(samples2[,2], rate = params2)

# Estimación de la densidad
dens1 <- kde(as.matrix(samples1)) #  dependencia positiva
dens2 <- kde(as.matrix(samples2)) #  dependencia negativa

# Convertir la estimación de densidad a un formato adecuado
x1 <- dens1$eval.points[[1]] 
y1 <- dens1$eval.points[[2]]
z1 <- matrix(dens1$estimate, nrow = length(x1), ncol = length(y1))

x2 <- dens2$eval.points[[1]] 
y2 <- dens2$eval.points[[2]]
z2 <- matrix(dens2$estimate, nrow = length(x2), ncol = length(y2))

# Crear gráficos interactivos con plotly
p1 <- plot_ly(x = ~x1, y = ~y1, z = ~z1) %>% add_surface()
p2 <- plot_ly(x = ~x2, y = ~y2, z = ~z2) %>% add_surface()

2.1.4.2.1 Caso Positivo

Código

p1

2.1.4.2.2 Caso Negativo

Código

p2

2.1.5 Arquimediana Frank

\[C(u,v)=\phi^{-1}(\phi(u)+\phi(v))\]

donde \(\phi(u)=-\ln\frac{e^{-\theta u}-1}{e^{-\theta }-1}\)

2.1.5.1 Función de distribución conjunta

\[H(x, y) = -\frac{1}{\theta} \ln \left( 1 + \frac{(e^{-\theta (1 - e^{-\lambda x})} - 1)(e^{-\theta (1 - e^{-2\lambda y})} - 1)}{e^{-\theta} - 1} \right)\]

2.1.5.2 Función de densidad conjunta

Código

theta1 = 10 #  dependencia positiva
theta2 = -15 #  dependencia negativa

# Definir cópulas de Clayton
myCopula1 <- frankCopula(theta1, dim = 2)
myCopula2 <- frankCopula(theta2, dim = 2)

# Número de muestras
n <- 10000

# Generar muestras de la cópula
set.seed(123)
samples1 <- rCopula(n, myCopula1) #  dependencia positiva
samples2 <- rCopula(n, myCopula2) #  dependencia negativa

# Convertir las muestras a distribuciones exponenciales
samples1[,1] <- qexp(samples1[,1], rate = params1)
samples1[,2] <- qexp(samples1[,2], rate = params2)
samples2[,1] <- qexp(samples2[,1], rate = params1)
samples2[,2] <- qexp(samples2[,2], rate = params2)

# Estimación de la densidad
dens1 <- kde(as.matrix(samples1)) #  dependencia positiva
dens2 <- kde(as.matrix(samples2)) #  dependencia negativa

# Convertir la estimación de densidad a un formato adecuado
x1 <- dens1$eval.points[[1]] 
y1 <- dens1$eval.points[[2]]
z1 <- matrix(dens1$estimate, nrow = length(x1), ncol = length(y1))

x2 <- dens2$eval.points[[1]] 
y2 <- dens2$eval.points[[2]]
z2 <- matrix(dens2$estimate, nrow = length(x2), ncol = length(y2))

# Crear gráficos interactivos con plotly
p1 <- plot_ly(x = ~x1, y = ~y1, z = ~z1) %>% add_surface()
p2 <- plot_ly(x = ~x2, y = ~y2, z = ~z2) %>% add_surface()

2.1.5.2.1 Caso Positivo

Código

p1

2.1.5.2.2 Caso Negativo

Código

p2

2.1.6 Arquimediana Ali-Mikhail-Haq (AMH)

\[C(u,v)=\phi^{-1}(\phi(u)+\phi(v))\]

donde \(\phi(u)=\ln\frac{1-\theta(1-u)}{u}\)

2.1.6.1 Función de distribución conjunta

\[H(x, y) =\frac{(1 - e^{-\lambda x})(1 - e^{-2\lambda y})}{1 - \theta (e^{-\lambda x})((e^{-2\lambda y}))} \]

2.1.6.2 Función de densidad conjunta

Código

theta1 = 1 #  dependencia positiva
theta2 = -1 #  dependencia negativa

# Definir cópulas de Clayton
myCopula1 <- amhCopula(theta1, dim = 2)
myCopula2 <- amhCopula(theta2, dim = 2)

# Número de muestras
n <- 10000

# Generar muestras de la cópula
set.seed(123)
samples1 <-  rCopula(n, myCopula1) #  dependencia positiva
samples2 <-  rCopula(n, myCopula2) #  dependencia negativa

# Convertir las muestras a distribuciones exponenciales
samples1[,1] <- qexp(samples1[,1], rate = params1)
samples1[,2] <- qexp(samples1[,2], rate = params2)
samples2[,1] <- qexp(samples2[,1], rate = params1)
samples2[,2] <- qexp(samples2[,2], rate = params2)

# Estimación de la densidad
dens1 <- kde(as.matrix(samples1)) #  dependencia positiva
dens2 <- kde(as.matrix(samples2)) #  dependencia negativa

# Convertir la estimación de densidad a un formato adecuado
x1 <- dens1$eval.points[[1]] 
y1 <- dens1$eval.points[[2]]
z1 <- matrix(dens1$estimate, nrow = length(x1), ncol = length(y1))

x2 <- dens2$eval.points[[1]] 
y2 <- dens2$eval.points[[2]]
z2 <- matrix(dens2$estimate, nrow = length(x2), ncol = length(y2))

# Crear gráficos interactivos con plotly
p1 <- plot_ly(x = ~x1, y = ~y1, z = ~z1) %>% add_surface()
p2 <- plot_ly(x = ~x2, y = ~y2, z = ~z2) %>% add_surface()

2.1.6.2.1 Caso Positivo

Código

p1

2.1.6.2.2 Caso Negativo

Código

p2

2.2 Caso Pareto

Tenemos las siguientes variables aleatorias:

\[X\sim \textbf{Pareto}(\alpha, 2\nu)\]

\[Y\sim \textbf{Pareto}(\alpha, \nu)\]

donde \(F_{x}(x)=1-\left(\frac{2\nu}{x+2\nu}\right)^{\alpha}\) y \(F_{y}(y)=1-\left(\frac{\nu}{y+\nu}\right)^{\alpha}\)

2.2.1 Cópula Gaussiana

\[C(u,v)=\Phi_{2}(\Phi^{-1}(u),\Phi^{-1}(v))\]

2.2.2 Cópula Marshall-Olkin

\[C(u,v)=(1-u)(1-v)\textbf{mín}\{(1-u)^{-\alpha_{1}},(1-u)^{-\alpha_{1}}\}\]

2.2.3 Arquimediana 1

\[C(u,v)=\phi^{-1}(\phi(u)+\phi(v))\]

donde \(\phi(u)=(-\ln u)^{\theta}\)

2.2.4 Arquimediana 2

\[C(u,v)=\phi^{-1}(\phi(u)+\phi(v))\]

donde \(\phi(u)=\ln(1-\theta\ln(u))\)

2.2.5 Arquimediana 3

\[C(u,v)=\phi^{-1}(\phi(u)+\phi(v))\]

donde \(\phi(u)=(1-\ln(u))^{\theta}-1\)